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We show theoretically that an open-dissipative polariton condensate confined within a trapping 
potential and driven by an incoherent pumping scheme gives rise to bistability between odd and even 
modes of the potential. Switching from one state to the other can be controlled via incoherent pulsing 
which becomes an important step towards construction of low-powered opto-electronic devices. The 
origin of the effect comes from modulational instability between odd and even states of the trapping 
potential governed by the nonlinear polariton-polariton interactions. 


PACS numbers: 71.36.+C, 42.65.Pc, 42.55.Sa 

I. INTRODUCTION 

Switching waves arising in bistable light-matter systems 
are of fundamental importance for future opto-electronic 
devices that intertwine both the characteristics of elec¬ 
trons and photons. Such waves can be thought of as mov¬ 
ing domain walls (wave crests) where one state succumbs 
to its bistable counterpart. In general, light-matter sys¬ 
tems owe their bistability to their optical non-linear na¬ 
ture. An example of such a system is the optical mi¬ 
crocavity 1 where strong coupling between excitons living 
in embedded quantum wells, and cavity photons gives 
rise to the light exciton-polariton quasiparticle (or sim¬ 
ply, polariton) 2 . In the past decade, a great amount of 
research has been dedicated to these light-matter bosons. 
Among the most interesting properties of the polariton 
is its ability to condense into a macroscopically occupied 
coherent state 3 , essentially a Bose-Einstein condensate, 
at temperatures much higher than conventional atomic 
condensates 4 due to its light effective mass. 

It has been shown that exciton-polariton systems give 
rise to bistability under resonant pumping schemes 5-7 
with applications in photonic logic gates 8-10 , switches 
and memory elements 11-13 , control of polariton superflu¬ 
ids 14 , and optical diodes 15 . Under resonant excitation a 
number of works focused on the parametric instability be¬ 
tween polaritons 16 20 , where polaritons in one particular 
energy-momentum state scatter in pairs to ” signal” and 
” idler” states conserving energy and momentum. Theo¬ 
retically, a bistability of the system was shown to be asso¬ 
ciated with the process 7 , corresponding to the stability 
of states under the same conditions where the scatter¬ 
ing could be activated or not. Bistability resulting from 
parametric interaction was further shown to allow for the 
formation of solitons 21 23 . 

To lift the requirement of a resonant laser, a few works 
have considered mechanisms to generate bistability under 
non-resonant or incoherent excitation, but typically re¬ 
quired specially designed structures 24 . Theoretical work 
showed that parametric instability and an associated 
bistability could be arranged in sub-wavelength grating 


microcavities 25 . Recently, the existance of spin-based 
bistability induced by different XY -polarization lifetimes 
was shown experimentally 26 . It has also been shown that 
modulational instabilities can lead to new states that co¬ 
exist with stable states 27 . 

In this work we study a mechanism of bistability, 
appearing when polaritons are incoherently excited in 
the presence of a trapping potential. Different modes 
of the trapping potential can be preferentially excited 
via the spatial patterning of the incoherent excitation. 
Non-linear coupling between the modes generated by 
polariton-polariton interactions then leads to a modu¬ 
lational instability. Remarkably, a bistability is associ¬ 
ated with this process and switching between different 
bistable states can be achieved with incoherent pulses, 
even when only the first two eigenmodes are considered. 
In extended 2D systems, switching waves that determine 
the boundary between spatially separated domains can 
be excited. 

The article is structured as follows: In Sec. II we in¬ 
troduce the model describing the coherent polariton con¬ 
densate confined within a trapping potential. In Sec. Ill 
we present results where we limit our analysis to only the 
first two eigenstates of the harmonic potential. In Sec. IV 
we present results considering only the first three eigen¬ 
states of the harmonic potential. In Sec. V we present 
results considering the full theoretical model containing 
all eigenstates of the harmonic potential. Here we demon¬ 
strate a switching wave signal traveling in a 2D planar 
system. In Sec. VI we demonstrate analogous bistability 
in another type of trapping geometry, the infinite quan¬ 
tum well. 


II. THEORETICAL MODEL 

We consider a microcavity with an embedded 2D quan¬ 
tum well subjected to a trapping potential 28,29 . In an 
open-dissipative system where the incoherent pumping 
is balanced against the polariton lifetime and an exciton 
reservoir decay mechanism 30 , the true solutions of the in¬ 
teracting condensate confined in a trapping potential are 
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Figure 1. (Color online) Example schematic showing the case 
of a harmonic potential, V(x) = muj 2 x 2 /2 and its first three 
states. The incoherent pump P(x) (orange) is chosen to have 
maximum overlap with one state but will also have some fi¬ 
nite overlap with higher states thus creating a condensate of 
coupled modes. 


where Pq denotes the intensity of the pump. While the 
pump is chosen so as to have a maximum overlap for the 
A\ mode making it condense first, it understandably has 
non-zero overlap with other modes (see Fig. 1). This re¬ 
sults in an interplay between the gain and the decay of 
all the modes. Substituting Eqs. 2-3 into Eq. 1 and inte¬ 
grating over the n-th state, i/j n , we come to the dynamical 
equation for the condensate modes, 

dA n (_ . v—> 

= V ^ n 2t ) 3 ^ > MnllmA-m . . . 

' ' m 

+ (a- iR ) ]T M njkl A*A k At. (4) 

jkl 

Here are transition matrix elements that obey, 

/ oo 

'ipi'ipjipki’i dx. (5) 

-OO 

III. TWO-MODE MODEL 


in fact non-trivial 31 . The trapping potential eigensolu- 
tions (modes) in the ID Schrodinger equation are termed 
'ipn(x) with energies E n and form a complete basis such 
that we can always write our condensate order param¬ 
eter in this basis. Assuming that the eigenenergies E n 
are different, then in order for the condensate to become 
stationary, populations in different modes are blueshifted 
such to become degenerate, allowing polaritons freely to 
’flow’ from one mode to another. 

To account for the aforementioned effects, the spatial 
dynamics of interacting polaritons can be described by 
the complex Gross-Pitaevskii (cGP) equation 32 , 
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h 2 d 2 
2m dx 2 
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(i) 


Here T is the polariton condensate order parameter, m is 
the polariton effective mass, a is the interaction strength, 
r is the polariton lifetime, R is the condensate saturation 
(also known as reservoir depletion rate), P(x) is the spa¬ 
tial profile of the incoherent pump intensity, and V(x) is 
the trapping potential. 

The condensate order parameter can be written, 


(x, t) = ^2 A n (t)lpn ( X ). (2) 

n=0 

The time evolution of the cGP-equation is governed by 
the modes amplitudes A n (t). In order to favor the gen¬ 
eration of polaritons in the first excited state the pump 
takes the following profile, 

P(x) = Po|V>i(s)| 2 (3) 


In order to test our theory for the most simple case, we 
investigate whether bistability can exist considering only 
scattering between the first two states of the ID harmonic 
oscillator, ipo, and t/q, assuming that modulational insta¬ 
bility takes place. Modulational instability was also pre¬ 
dicted for non-resonantly excited polariton condensates 
where the condensate has a strong back-action effect on 
an incoherent reservoir of hot exciton states 33,34 . We 
are not in this regime as the complex Gross-Pitaevskii 
approach 32 , corresponding to Eq. 1, applies to the case 
where the reservoir has a faster timescale than polariton 
dynamics and can be treated as independent. 

The condensate order parameter can be written in the 
limited basis of the first two harmonic oscillator states, 

^(x,t) = A 0 (t)ip o(x) + Ai(i)^>i(a:). (6) 

A complete analysis is given in appendix A. We begin 
by adiabatically switching on the pump. A condensa¬ 
tion threshold is reached for the A\ mode which, in 
the low density regime, satisfies Pq = h/2rMiin. At 
higher pumping powers we find that indeed there exists 
a bistable regime (see Fig. 2) where switching from one 
solution to the other can take place by pulsing the system 
incoherently. The blue and red solid (dashed) lines show 
results of slowly increasing (decreasing) the pump. The 
green and black dots are semi-analytical results given by 
Eq. A15 in appendix A. One can see that the central 
branches of the hysteresis curve are resolved from the 
Eq. A15 whereas numerically time resolved results only 
show the bottom and the top branches. 

In the case of a noninteracting condensate where a = 0 
and setting r 4 oo we find stationary solutions corre¬ 
sponding to the case where only one mode stays popu¬ 
lated (see appendix A). Numerically we observe no so¬ 
lutions where both modes are populated at the same 
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Figure 2. (Color online) Whole (dashed) blue and red lines 
show the time resolved density from Eqs. A1-A2 when the 
pump is slowly switched on (off). Calculated solutions from 
the analysis leading up to Eq. A15 (green and black dots). 
Parameters: m — 3 • 10 _5 mo, hco — 51.5 /xeV, a — 2.4 /xeV 
/im, r = 2 ps and R — 0.3a. 

time. Indeed, from Eqs. A1-A2 it is not directly evi¬ 
dent whether such a solution exists or not. Thus, in 
theory, it remains unproven whether bistability exists in 
the absence of interactions. However, in Sec. V[a] we 
will demonstrate that the full noninteracting cGP-model 
yields population higher modes. 



Figure 3. (Color online) Results of the three-mode model 
under adiabatic incoherent pumping, (a) Population of the 
modes shown when the pump is slowly switched on (solid 
lines) and then slowly turned off (dashed lines). At a certain 
pump threshold there is a transition from A\ into a conden¬ 
sate composed of all three modes. At higher pump values only 
Ao and A 2 stay populated. One can see that a bistable inter¬ 
val exists at 2.1 meV /im < Po <3.1 meV /un. (b) A close up 
of the transition from frame (a). Parameters: m = 3T0 -5 rao, 
huj — 34 /xeV, a — 2.4 jueV jum, r = 15 ps and R — a. 


IV. THREE-MODE MODEL 

We now extend our model to the first three eigenmodes 
of the HO where the order parameter is written, 

^(x,t) = A 0 (t)ip 0 (x) + Ai(t)ipi(x) + A 2 {t)ij) 2 {x). (7) 

Solving the three mode model numerically (see Eqs. BI¬ 
BS) we find that bistable regions can indeed exist through 
a wide range of oscillator frequencies u and pump inten¬ 
sities Po (see Figs. 3 and 4). Results displayed in Fig. 3 
show the change in population of the three modes as the 
non-resonant pump is slowly switched on (whole lines) 
and then slowly switched off (dashed lines). First, a con¬ 
densation threshold is reached for the A\ mode which in 
the low density regime satisfies Po = h/2rMnn ~ 0.65 
meV /im. At higher pumping values there is an abrupt 
drop in the strength of the A\ mode as scattering and 
gain of the Ao and A 2 modes overcome their decay. Just 
as with the two-mode model, the three modes are degen¬ 
erate and co-exist at the same energy at all times when 
they are populated. As the pump is increased further, 
the A 1 mode is completely quenched and only mode Ao 
and A 2 stay populated. 

When the pump is slowly decreased (dashed lines in 
Fig. 3), one can clearly see that modes A 0 and A 2 re¬ 
main supported over an interval of pump values where 
previously only the A\ mode was supported. A phase di¬ 
agram showing the bistable region across different values 
of oscillator energies Jtlj and pump intensities Po is dis¬ 
played in Fig. 4. The black dashed lines indicate the two 
different bistable areas. The larger one, (i), shows where 
only A\ exists but can be excited to a solution where Ao 


and A 2 only exist. The smaller one, (ii), shows where all 
modes are populated but can also be excited to an Ao 
and A 2 only solution. 

Analogous to our results in Sec. Ill, numerically we ob¬ 
serve no population in the Ao and A 2 modes when inter¬ 
actions are absent. 
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Figure 4. (Color online) Population of the second mode A\ 
shown for varying oscillator frequency and pump intensity. 
Dashed lines enclose the two different bistable regions corre¬ 
sponding to those demonstrated in Fig. 3. Same parameters 
are used as in Fig. 3. 


V. FULL CGP-FORMALISM 

So far simplified models have been considered using only 
the first two- and three eigenmodes of the HO. These 
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Figure 5. (Color online) (a) Populations of the first eight 
harmonic modes using the full cGP formalism. A bistable in¬ 
terval is observed at 9.8 meV /im < Po < 13.9 meV /im when 
pump is increased (solid lines) and then decreased (dashed 
lines) associated with odd modes transitioning into even ones 
and back, (b) Condensate energy during the forward tran¬ 
sition for different potential strengths, (c) Increasing the 
pump strength reveals that the transition from odd to even 
modes takes place periodically with increasing pump power 
and does not affect the total integrated density Np (black 
line). Note that frames (b-c) are only shown for increasing 
pump strength. Parameters were set to: m — 3 • 10 _5 mo, 
r = 2 ps, a = 2.4 /ieV /im and R = 0.3o. In frames (a,c): 
ftw = 30 /reV. 


This is analogous to the results displayed in Fig. 3 where 
A\ succumbs to Ao and A 2 . Slowly decreasing the pump 
(dashed lines in Fig. 5 [a]) reveals a bistable region be¬ 
fore the condensate transitions back to a superposition of 
odd modes. Both solutions of the condensate are found 
to be stable within 9.8 meV fi m < Po ^ 13.9 meV /im. 
From Fig. 5 [a] it can be seen that the population strength 
of different modes does not follow a strict hierarchy, an 
example is the crossover of |As| 2 and \A?\ 2 . The en¬ 
ergy of the condensate reveals a sudden blueshift when 
the transition takes place (see Fig. 5 [b]) which increases 
with the oscillator energy hu. Interesting enough, the ob¬ 
served transition takes place periodically with increasing 
pump intensity (see Fig. 5 [c]), and using pump spatial 
profiles corresponding to higher harmonic states also in¬ 
duces a transition from and odd to an even condensate. 
In Fig. 5[c] we see that the total integrated density Np 
(black line) is unaffected by these transitions. We fur¬ 
thermore observe phase-locking taking place between the 
two condensate solutions as the condensate is transition¬ 
ing from one to another, corresponding to coherent trans¬ 
port of polaritons between odd and even states. Such 
phase synchronization phenomenon has been previously 
reported in spinor polariton condensates due to para¬ 
metric scattering 35 . In the noninteracting case (a = 0) 
the transition is no longer observed for increasing pump 
power and the polariton population remains only in the 
odd modes, corresponding to the nontrivial stationary 
state of the pumped noninteracting condensate. This is 
in analogy with results from Secs. Ill-IV where no popu¬ 
lation could be observed in the neighboring even modes, 
Aq and A 2 . 


have demonstrated two artifacts of the trapped incoher¬ 
ently generated polariton condensate. Namely, a bistable 
interval between the mode populations (Sec. III-IV), and 
a sudden transition from an odd solution into an even 
one (Sec. IV). 

Realistic polariton systems however would contain many 
spatial modes, corresponding to a full model of the polari¬ 
ton spatial dynamics according to Eq. 1. Here, we solve 
the Eq. 1 directly by propagating it stepwise in time and 
demonstrate that both of these artifacts from previous 
sections are still clearly present. 


A. ID HO System 

Results for a ID system are shown in Fig. 5, calcu¬ 
lated for typical parameters of a GaAs based system: 
m = 3 • 10 _5 mo, where mo is the free electron mass, 
a = 2.4 /ieV /im, r = 2 ps and R = 0.3<x As the 
pump intensity is slowly increased we observe population 
growth only in the odd modes with the fastest growth in 
A\. As the pump intensity increases further a sudden 
transition takes place where the population of the odd 
modes (Ai included) suddenly shifts to the even modes. 


B. 2D HO System 

Applying the ID system to a 2D one, the polariton con¬ 
densate can be confined along a parabolic guide. By 
slowly switching on the non-resonant background pump 
to the bistable regime, a nice steady state solution, su¬ 
perposed of odd modes, is formed along the guide. The 
condensate is then pulsed incoherently at one end result¬ 
ing in a switching wave which travels along the guide (this 
can also be achieved by coherent injections) transitioning 
odd harmonic modes into even modes. In Fig. 6(a) we see 
a signal traveling along a straight guide, set to fruo = 60 
/ieV, after a pulsed incoherent injection on the left end. 
At first the signal moves slowly, then it speeds up as the 
condensate start to switch rapidly into its bistable coun¬ 
terpart. From numerical calculations the velocity of the 
signal is determined to be around ~ 1.2 fi m/ps. 

The temporal width of the pulse is in the order of a few 
picoseconds. Fig. 6(b) shows the corresponding change 
in iF-space occupancy, extracted from x = 0 in Fig. 6(a). 
In Fig. 6(c) we demonstrate the same type of a switching 
wave but now in a guide with huo = 30 /ieV. Here the sig¬ 
nal travels successfully around a curved guide completing 
a 90 degree turn. Such switching waves can be realized 
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Figure 6. (Color online) (a) Switching wave demonstrated in 
a 2D system with a parabolic trap of strength ftw — 60 /xeV 
along the x-axis. An incoherent pulse is activated causing a 
sudden transition from a condensate composed of odd har¬ 
monic states into even states in a form of a switching wave, 
(b) K -space dynamics showing change of topology when the 
transition occurs, (c) A switching wave traveling along curved 
potential geometry where fouj = 30 /xeV. Other parameters 
were set to same values used in Fig. 5. 


over a large range of potential frequencies. In fact, in¬ 
creasing the trap frequency u reveals that the blueshift 
associated with the switching increases (Fig. 5 [b]) and 
consequently increases the speed of the signal This opens 
up the possibility of controlling the signal speed using 
different potential strengths. 
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Figure 7. (Color online) Bistable area demonstrated in the 
infinite quantum well of width L = 60 /xm, the result is anal¬ 
ogous to the bistable area displayed in Fig. 5. Parameters: 
a — 2.4 /xeV pm, R — 0.3a:, r = 2 ps, m = 3 • 10 _5 mo. 


VI. INFINITE QUANTUM WELL 

Another example of a well known trapping geometry 
is the infinite quantum well, a well studied system with 
known eigensolutions to the Schrodinger equation. The 
main difference here from the HO is that the energy levels 
are no longer equidistant but grow quadratically with the 
mode quantum number n. Applying the same method as 
in Sec. V[a] we see that a bistable interval also exists (see 
Fig. 7). Note the indexing of the modes in Fig. 7 follows 
common literature on the infinite quantum well, where 
the quantum number starts at n = 1 whereas in the HO 
it starts at n = 0. 


VII. CONCLUSIONS 

We have shown that an open-dissipative polariton con¬ 
densate, supported by incoherent pumping, described by 
a Gross-Pitaevskii type equation confined in a trapping 
potential can possess bistable regions. The bistability 
is inherited from firstly scattering between degenerate 
condensate modes, and secondly gain-decay mechanisms 
resulting from the different pump overlap for different 
harmonic modes. Results show that the minimum re¬ 
quirement for such bistability to take place is a system 
containing only two separate energy levels. 

For higher pumping powers the condensate undergoes 
a series of oscillations in parity between highly occupied 
odd and even modes. We can switch from one solution 
of the condensate to its bistable counterpart in 2D har¬ 
monic potential guides, carrying information at veloci¬ 
ties around ^ 1.2 /xm/ps, and managing to travel along 
curved guides. 
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Appendix A: Two-Mode Model 

Starting from Eq. 4, we come to the following two coupled equations for the first two eigenstates of the harmonic 
oscillator (neglecting higher modes) described by the cGP equation. 

ih = (e 0 — —^ Aq + iPoMoonAo + (a — iR) [M 0 ooo|Ao| 2 A 0 + Mxm (AqA 2 + 2A 0 |Ai| 2 )] (Al) 

**£ = (Eh. iPoMnnAj. + (a- iR) [M lul \A 1 \ 2 A l + M 00 ii (2|Ao| 2 A! + A 2 A \)] . (A2) 

Let’s keep in mind that the matrix elements can be written in a more simpler form, 

I mu /3 3/3 

Afoooo = J ^oon = —, Afnn = —. 


1. Noninteracting Case 


Though the coupled dynamical equations offer four independent equations to solve the problem exactly for any finite 
population in both modes, the nonlinear nature of the coupling does not reveal a clear solution depending on the 
physical parameters of the problem. It remains then unproven whether for some set of (uj, Po, R ) one can have solution 
with finite populations in both modes in the absence of interactions. 

However, one can determine two trivial stationary solutions: |Ao| 2 = 0, |Ai| 2 = Po/P, and |Ao| 2 = Po/2P, \Ai\ 2 = 0. 
In princible, one can then have bistability if both solutions are stable. To check this, we perform stability analysis 
analogous to Ref. 36 on the two solutions. We use the following ansatz for the former solution, 

Mt) = e~ ilxt (A° 0 + SA 0 (t )) , A x (t) = e"** (0 + SA^t)) 


Let’s plug into Eqs. A1-A2 and neglect higher order terms of <SAo,i. The steady state will satisfy 
h/i = Eq. 


. k d{5A o) , n /y 4 . no (SAo)* 

in —--= —iPq—oAq — iPqp -. 

dt 2 P 2 ' 


Po/2R and 


(A3) 


ih 


djSAj) 

dt 


(hio + iP 0 Aj 


SAi — iPo 


P(SAiF 
2 2 


(A4) 


The fluctuations can be described with <5A(yi (t) = Mqj e l£l + 1 e l£ 1 where he is the fluctuation energy. We then 

come to an eigenvalue problem = chT where <5T := (uq,vq, ui,v\) t . The matrix is, 


(-iPoP/2 -iPoP/2 0 

-iP 0 p/2 —iPo/3/2 0 

0 0 hui + iPo/3/4: 

\ 0 0 -iPoP/4 


0 

0 


-iPoP/* 

—hu + zPq/3/4, 


(A5) 


and for the case of A 0 (£) = e lflt (0 + SA 0 (t)) and Ai(t) = e (AJ + SAi(t)) we have, 


/-hu - iPo/3/2 -iPo/3/4 0 ° \ 

-iPo/3/4: hu-iP 0 f3/2 0 0 

0 0 -i3P 0 f3/4 -i3P 0 (3/8 

\ 0 0 -i3P 0 f3/8 —i3Pof3/4:J 


(A6) 


The spectrum will then consist of four eigenvalus Si corresponding to the normal modes of 5Ao,i. Here the stability 
is determined by the imaginary part of the eigenvalues. If all imaginary parts are negative the fluctuations decay 
exponentially. If however one or more eigenvalue is positive the the fluctuations grow exponentially. In Fig. 8 we plot 
the imaginary parts of eigenvalues for Eqs. A5-A6 for different values of Po and uj. We find that positive imaginary 
parts of the eigenvalues exist when Aq is populated underlining that it is unstable whereas when we have only 
population in A\ it is stable. 
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Figure 8. (Color online) (a-b) Imaginary parts of the eigenvalues given by Eq. A5-A6 respectively for varying oscillator frequency 
uj and Pq = 1 meV /im. (c-d) Same for ftw — 30 /xeV and varying pump intensity Pq. For all plots: R = 0.72 /xeV /xm. 


2. Interacting Case 

In order to find an analytical expression for the two modes when they are both populated and interactions are 
nonzero, we make an attempt at solving Eqs. A1 and A2 for some arbitrary complex solution Ao = (a + and 

A\ — ce _2Axt where we have assumed that the two modes are degenerate with condensate energy hp and the global 
phase is invariant (only relative phase is important). We then arrive at, 


, huo ih 

hp = — - — + (a-iR) 


M 0 ooo (& 2 + b 2 ) + Mqoiic 2 


2 iab 


b 2 


+ 2Moonc 2 


hfi — ——-——h (o' — iR) Mime 2 + Moon (a 2 — b 2 + 2iab) + 2Moon(tt 2 + b 2 ) 

2 2r . 


+ iPoMoon 
+ iPoMmi 


(A7) 

(A8) 


From the first equation: 


^ /9 9 . „ 9 a 2 — b 2 „ 9 _ ^ afrc 2 

0 = < —- h/i + a/3(a 2 + b 2 ) + a/3c 2 o/r 9 | + <a/3c 2 — i?/3- 


2(a 2 + 6 2 ) 


a 2 + b 2 


. r n 

— i. < — 


ab 


'{~2~ P +a ^^i 


2 7.2 

■ Rf3(a 2 + 6 2 ) + + We? ~ p o 

2(a 2 + b z ) 


Re 

-Ml 


2 J 


Im 


(A9) 

(A10) 


From the second equation: 


f 3 Huj 

0 — -hfi + a 4 


+ a^(a 2 — b 2 ) + a/3(a 2 + b 2 ) + i?/3a&]- 


-{A 


3/3 


a/3ab + R^~c z + i7/3 


) Re 

+ RP(a? + b 2 )-pM\ 

4 J Im 


(All) 

(A12) 


The real and imaginary parts must be zero. Let’s take the difference of the real parts and imaginary parts respectively 
to get two separate equations for c 2 and to get rid of fi and r for the time being. 


2 2 a(a 4 — b 4 ) + 4i7a6(a 2 + b 2 ) + 4(/kj//3)(a 2 + b 2 ) 


c 2 = 


2a(a 2 — b 2 ) + a(a 2 + b 2 ) — 4Rab 
2 R(a 4 — b ‘ 4 ) — 4 aab(a 2 + b 2 ) — Po(a 2 + b 2 ) 

2 R(a 2 - b 2 ) + R(a 2 + b 2 ) + 4aab E 


(A13) 


(A14) 























Equating the two, we arrive at a cubic equation for b, 


0 = -6 3 [3 a(a 2 + R 2 )] - b 2 


hoj P 0 

T fl + T“ 


+ b 


Ahuj 


5 a 6 (a 2 + it ) H-— aa — PqRcl 


+ 


3hu) 

~T 


t~) 2 | 3Pp 2 
Ra H-— aa 


(A15) 


Thus, the bistable area corresponds to an interval where three real roots exist to this equation for a real valued c and 
satisfy the constraint that Eqs. A9-A12 should be zero. This allows us to solve the hysteresis branches of bistable 
areas as shown in Fig. 2 (green and black dots). 

Appendix B: Three-Mode Model 

Working with only the first three eigenstates of the HO. One arrives at three coupled equations, 
iti dt = ~ 2r / ^o(AfoonAo + A/ 0112 A 2 ) + (a — iR)\M qqqq\Aq\ 2 Aq + A^fo 2221 ^- 21 2 ^-2 • • • 


+ A /0002 


2|A 0 | A 2 + AqA% 


+ Afoon 


AqA( + 1A 0 |Ai| 


+ AfoQ22 


AqAz + 2A 0 |A 2 |" 


+ A /0112 


2|Ai|A4 2 + AfA* 


(Bl) 
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